

%% load estimation results
load('parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
clearvars -except thetaStar K mpow N reb X Y YY Z specif


%% load everything of interest
load('sim_eq_round2.mat')
[DMES,intervEq] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);


%% standard errors via Hessian
InvSigmaH = HessPhat(thetaStar,1,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow);
seH = sqrt(diag(inv(InvSigmaH)));
resultsH = printResults([thetaStar,seH,2*normcdf(abs(thetaStar./seH),'upper')],reb,mpow,specif);
disp(' ')
disp('for Tables 1, 2, and D1:')
disp(resultsH)


%% standard errors via parametric bootstrap

% "long" replication
% load('model_fit_replication.mat','Yh')
% load('inference_replication.mat','srngBoot','B')
% thetaBoot = zeros(length(thetaStar),B);
% exitBoot = ones(1,B);
% timeBoot = 0;
% parfor b = 1:B
%     tic
%     rng(srngBoot(b))
%     Yb = arrayfun(@(n)datasample(1:size(YY,1),1,'Weights',Yh(:,n)),1:N)'; %#ok<PFBNS>
%     optionsb=optimset('Display','off','HessFcn',@(x,lambda)HessPhat(x,lambda,Yb,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow));
%     theta0=thetaStar+[zeros(length(thetaStar),1),randn(length(thetaStar),2).*abs(thetaStar)/2]; % starting values
%     l0=1:3; exit0=1:3;
%     for sv = 1:3
%         l0(sv) = Phat(theta0(:,sv),Yb,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow);
%         if ~isfinite(l0(sv))
%             theta0(:,sv) = NaN(length(thetaStar),1);
%             l0(sv) = NaN;
%             exit0(sv) = NaN;
%         else
%             [theta0(:,sv),l0(sv),exit0(sv)]=knitromatlab(@(x)Phat(x,Yb,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow),...
%                 theta0(:,sv),[],[],[],[],[],[],[],[],optionsb,'knitro1.opt');
%         end
%     end
%     [~,sv] = min(l0);
%     thetaBoot(:,b) = theta0(:,sv);
%     exitBoot(:,b) = exit0(sv);
%     timeBoot=timeBoot+toc;
% end
% seBoot1 = std(thetaBoot,0,2);
% load('inference_replication.mat')
% disp(' ') 
% disp('verify bootstrap results:') 
% disp(['max. rel. diff. between std. errors = ',num2str(max(abs(seBoot-seBoot1)./seBoot))])
% resultsBoot = printResults([thetaStar,seBoot,2*normcdf(abs(thetaStar./seBoot),'upper')],reb,mpow,specif);

% "short" replication
load('inference_replication.mat')
disp(' ')
disp('for Table D1 (bootstrapped standard errors):')
disp(resultsBoot)


%%
clearvars -except thetaStar InvSigmaH seH resultsH srngBoot B ...
    thetaBoot thetaBoot0 exitBoot timeBoot seBoot resultsBoot
save('inference.mat')



